The hydration number of Na + in liquid water 
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An 'ab initio' molecular dynamics simulation of a Na + ion in aqueous solution is presented and 
discussed. The calculation treats a Na + ion and 32 water molecules with periodic boundary con- 
ditions on a cubic volume determined by an estimate of zero partial molar volume for this solute 
in water at normal density and at a temperature of 344±24 K. Analysis of the last half of the 
12 ps trajectory shows 4.6 water molecules occupying the inner hydration shell of the Na + ion on 
average, with 5 being the most probable occupancy. The self-diffusion coefficient observed for the 
Na + is l.OxlO" 5 cm 2 /s. The quasi-chemical theory of solutions provides the framework for two 
more calculations. First a complementary calculation, based on electronic structure results for ion- 
water clusters and a dielectric continuum model of outer sphere hydration contributions, predicts an 
average hydration shell occupancy of 4.0. This underestimate is attributed to the harmonic approx- 
imation for the clusters in conjunction with the approximate dielectric continuum model treatment 
of outer sphere contributions. Finally, a maximum entropy fitting of inner sphere occupancies that 
leads to insensitive composite free energy approximations suggests a value in the neighborhood of 
-68 kcal/mol for the hydration free energy of Na + (aq) under these conditions with no contribution 
supplied for packing or van der Waals interactions. 

keywords: ab initio molecular dynamics, dielectric continuum, electronic structure, hydration, 
information theory, quasi-chemical theory, sodium ion 



I. INTRODUCTION 

Solvation of simple ions in aqueous solution is not yet 
fully understood despite its fundamental importance to 
chemical and biological processes. For example, disagree- 
ment persists regarding the hydration number of the Na + 
ion in liquid water. A pertinent problem of current inter- 
est centers on the selectivity of biological ion channels; it 
seems clear that the selective transport of K + relative to 
Na + ions in potassium channels [Q, ||, || depends on de- 
tails of the ion hydration that might differ for K + relative 
to Na+. 

Experimental efforts to define the hydration structure 
of Na + (aq) using diffraction^, [j| and spectroscopic |6| 
methods produce a hydration number ranging between 
four and six[Q. Simulations have obtained a range of 
values, but most predict six water molecules in the in- 
ner hydration sphere of the Na+ ion|, § |(| |ll|, || 
Q [U| [0| g [19). An 'ab initio' molecular dynam- 
ics simulation produced five inner shell water molecules 
neighboring Na + (aq) Q ■ 

An important limitation of theoretical studies of ion 
hydration concerns the sufficiency of model force fields 
used in classical statistical mechanical calculations. In 
the most customary approaches, interatomic force fields 
used in theories or simulations derive from empirical 
fits of a parameterized model to a variety of exper- 
imental data. 'Ab initio' molecular dynamics avoids 
this intermediate modeling step by approximate solu- 
tion of the electronic Schroedinger equation for each con- 
figuration of the nuclei [pl| p2| . This technique thus 
goes significantly beyond conventional simulations re- 
garding the accuracy of the force fields. It also aug- 
ments theories built more directly on electronic struc- 
ture studies of ion- water complexes by adopting approx- 



imate descriptions of the solution environment of those 
complexes p3|, p4| , p5| , p(| |27| , |28|. 

Relative to conventional simulations, 'ab initio' molec- 
ular dynamics simulations also have some important lim- 
itations due to the high computational demand. Applica- 
tions of the method have been restricted to small systems 
simulated for short times. For example, an 'ab initio' 
molecular dynamics study pc[ of the Na + (aq) ion compa- 
rable to the present work obtained a thermal trajectory 
lasting 3 ps after minimal thermal aging. The present 
work, though still limited to relatively small systems, 
pushes such calculations to longer times that might per- 
mit more precise determination for Na + (aq) of primitive 
hydration properties. The analysis here utilizes the last 
half of a 12 ps thermal trajectory. The quasi-chemical 
theory^, |IJ |||, ||(| and separate electronic struc- 
ture calculations on Na(H20)„ + complexes assist in this 
analysis. 



II. METHODS 

The system consisted of one Na + ion and 32 water 
molecules in a cubic box with edge 9.86518 A and pe- 
riodic boundary conditions. The dimensions of the box 
correspond to a water density of 1 g/cm 3 and zero par- 
tial molar volume for the solute. Initial conditions were 
obtained as in an earlier 'ab initio' molecular dynamics 
simulation on Li + (aq) p3| . In that earlier work, an opti- 
mized structure for the inner sphere Li(Bi20)6 + complex 
was equilibrated with 26 water molecules under conven- 
tional simulation conditions for liquid water, utilizing a 
current model force field and assuming a partial molar 
volume of zero. In the present calculation, the same pre- 
equilibrated system was used as an initial configuration 
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for the 'ab initio' molecular dynamics except that an opti- 
mized structure for the inner sphere Na(H 2 0)g + complex 
replaced the hexa-coordinated Li + structure. Constant 
pressure or constant water activity simulations, defined 
by intensive rather than extensive variables, probably 
would produce a more useful characterization of the sol- 
vent thermodynamic state for these small systems, but 
those alternatives are currently impractical. 

Molecular dynamics calculations based upon a 
gradient-corrected electron density functional description 
of the electronic structure and interatomic forces were 
carried out on this Na + (aq) system utilizing the VASP 
program[^9). The ions were represented by ultrasoft 
pseudopotentials |30| and a kinetic energy cutoff of 31.5 
Rydberg limited the plane wave basis expansions of the 
valence electronic wave functions. The equations of mo- 
tion were integrated in time steps of 1 fs, which is small 
enough to sample the lowest vibrational frequency of wa- 
ter. A thermostat constrained the system temperature 
to 300 K during the first 4.3 ps of simulation time. After 
removing the thermostat, the temperature rose slightly 
and then leveled off by 6 ps to an average of 344 ± 24 K. 
During the simulation, the initial n=6 hydration struc- 
ture relaxed into n=4 and n=5 alternatives, such as those 
shown in Fig. [l| All analyses reported here rely on the 
trajectory generated subsequent to the 6 ps of aging with 
the system at a temperature elevated from room temper- 
ature. 



III. RESULTS 

The ion-oxygen radial distribution function is shown in 
Fig. |[ The first maximum occurs at a radius of 2.35 A 
from the Na + ion and the minimum at radius 3.12 A 
demarcates the boundary of the first and innermost hy- 
dration shell. An average of (n)=4.6 water molecules 
occupy the inner hydration shell. Fig. ^| tracks the in- 
stantaneous number of water oxygen atoms found within 
the first hydration shell of the Na + , defined by radius 
r<3.12 A for the upper panel. The fractions x^ and x§ 
of four-coordinate and five-coordinate hydration struc- 
tures, respectively, constitute ^4=40% and xs=56% of 
the last 6 ps of the simulation. Structures in which the 
Na + ion acquires six innershell water molecules occur 
with a 4% frequency, while structures with three and 
seven innershell water molecules occur less than 1% of 
the time. Analysis of the mean-square displacement of 
the Na + ion (Fig. |]) produces a self-diffusion constant of 
1.0 xl0~ 5 cm 2 /s, which agrees reasonably well with an 
experimental result of 1.33xl0~ 5 cm 2 /spl|. 

These results correspond coarsely with an 'ab initio' 
molecular dynamics calculation on this system carried- 
out independently ]20|| . The most probable inner shell 
occupancy found there was also five, but the probabil- 
ities of n=4 and n=6 were reversed from what we find 
here. This difference may be associated with the lower 
temperature used in Ref p0|. 



One motivation for this study arises from the quasi- 
chemical theory of solutions. According to this formu- 
lation, xq contributes a 'chemical' contribution to p e ^ a + , 
the excess chemical potential or absolute hydration free 
energy of the ion in liquid water p7), 



0fa 



= In xq — In 



(1) 



Here the inner shell is defined by specifying a function 
6 Na +j that is equal to one (1) when solvent molecule j is 
inside the defined inner shell and zero (0) otherwise; AU 
is the interaction energy of the solvent with the solute 
Na + that is treated as a test particle, /3 _1 =ksT, and 
the subscript zero associated with (. . .) indicates a test 
particle average [^7j . The second term on the right-hand 
side of Eq. (nh is the excess chemical potential of the 
solute lacking inner shell solvent molecules whereas the 
first term is the free energy of allowing solvent molecules 
to occupy the inner shell. The validity of Eq. (Q) has been 
established elsewhere p7[ . The second term on the right 
of Eq. ^ is the outer sphere contribution to the excess 
chemical potential in contrast to the first or chemical 
term. 

The utility of this quasi-chemical formulation is the 
suggestion p2[ of more detailed study of the x n , the frac- 
tions of n-coordinate hydration structures found in so- 
lution, on the basis of the equilibria forming inner shell 
complexes of different aggregation number: 



Na(H 2 O) m=0 + + nH 2 ^ Na(H 2 0) n - 
Utilizing the chemical equilibrium ratios 

„ PNa(H 2 0) n + 

Kn = n ' 

PH 2 PNa(H 2 O) m=0 + 

the normalized be expressed as 

_ K n pn 2 o n 
Xn ~ E K mP v 2 o m ■ 

m>0 



(2) 



(3) 



(4) 



The p a are the number densities and, in particular, pn 2 o 
is the molecule number density of liquid water. If the 
medium external to the clusters is neglected, the equi- 
librium ratios, denoted as K n (°\ can be obtained from 
electronic structure calculations on the complexes, as- 
suming for the thermal motion of the atoms the harmonic 
approximation evaluated at the calculated minimum en- 
ergy configuration. Finally utilization of a dielectric con- 
tinuum approximation for the outer sphere contributions 
to the chemical potential gives a natural, th oug h approx- 
imate, quasi- chemical model[^3|, ^J, |25|, g(| P4 |28||. 

For the present problem, the quasi-chemical approxi- 
mation was implemented following precisely the proce- 
dures of the earlier study of Li + (aq) , except that the 
sodium ion cavity radius for the dielectric model cal- 
culation was assigned as Rjv a +=3-l A, the distance of 
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the first minimum of the radial distribution function of 
Fig. ||. The temperature and density used were 344 K and 
1.0 g/cm 3 and the value of the bulk dielectric constant 
was 65.3@. 

Results of the calculations are summarized in Fig. [|. 
The electronic structure results are consonant with those 
found previously for the Li + ion. The n=4 inner sphere 
gas-phase complex has the lowest free energy. Although 
outer sphere placements are obtained for additional wa- 
ter molecules in the minimum energy structures of larger 
clusters, attention is, nevertheless, here restricted to in- 
ner sphere structures. The mean occupation number pre- 
dicted by this quasi-chemical model is (n) — 4.0; the com- 
puted absolute hydration free energy of the Na + ion un- 
der these conditions is -103 kcal/mol, not including any 
repulsive force (packing) contributions. An experimental 
value for Na + ion in liquid water at room temperature is 
-87 kcal/mol |5|. 

Because of the significance of x [Eq. , we fitted sev- 
eral model distributions {x n } based on different ideas to 
the 'ab initio' molecular dynamics results. The varying 
success of those models in inferring xq was enlightening. 
An instructive selection of those models is shown in Fig. o 
and we describe those results here. 

First, we note that though the preceeding quasi- 
chemical approximation does not agree closely with the 
'ab initio' molecular dynamics simulation, the popula- 
tions obtained from the quasi-chemical approximation, 
x n , can serve as a default model for a maximum entropy 
inference of x n f32]| . In this approach we model 

hiXj = IriXj - A -jXi - l)A 2 /2 - ... , (5) 

with Lagrange multipliers adjusted to conform to 
available moment information 




for j = 0, 1, 2, ... . In view of the limited data available, 
use of more than two moments produced operationally 
ill-posed fitting problems. 

One difficulty with this specific approach is that the 'ab 
initio' molecular dynamics produced xj>0 in contrast to 
the electronic structure methods that found no minimum 
energy hepta-coordinated inner-sphere clusters. Since 
the observed 27 is likely to be relatively less accurate 
and is furthest away from the desired n=0 element, we 
excluded n=7 configurations of the 'ab initio' molecu- 
lar dynamics, renormalized the probabilities x n , and re- 
calculated the moments. As the upper panel in Fig. ^ 
shows, this simple maximum entropy model is qualita- 
tively satisfactory although not quantitatively convinc- 
ing. The fitted model significantly disagrees with the 
observed 23. The chemical contribution suggested by 
Fig. ||is approximately -70 kcal/mol. Using the Born for- 
mula, -q 2 (l - l/e)/2R with R=3.12 A, to estimate the 
outer sphere contributions represented by the last term 
in Eq. [0, then the net absolute hydration free energy falls 



in the neighborhood of -115 kcal/mol. Since experimen- 
tal values for the absolute hydration free energy at room 
temperature center around -90 kcal/mol, this comparison 
shows that the present free energy results are not to be 
interpreted quantitatively, but rather as indicative of the 
present state of the theory. 

A second approach focused on testing a default model 
that supplies a nonzero £7; we used the Gibbs default 
model x n cx 1/n! that would give the correct answer for 
an ideal gas 'solvent.' This model has the additional 
and heuristic advantage of being significantly broader. 
Our experience has been that these maximum entropy 
fitting procedures work better when the default model is 
broader than the distribution sought. The results, illus- 
trated in the middle panel of Fig. show an improved 
fit. Here the chemical contribution to the free energy is 
-23 kcal/mol, yielding a net absolute hydration free en- 
ergy of -68 kcal/mol when the same Born formula is used 
to estimate the outer sphere contributions. 

A third fitting possibility was based on a suggestion 
from a previous 'ab initio' molecular dynamics calcula- 
tion on K + (aq): that the innermost four water molecules 
have a special status jj6|. In fact, the quasi-chemical ap- 
proximation above and the fitting of the upper panel of 
Fig. H suggests also that the x n results for n<4 and for 
n>5 display different behaviors. The radial distribution 
function of Fig. ||, somewhat better resolved than hereto- 
fore, is relevant to this issue and, in contrast, doesn't 
directly support a hypothesis of two populations of wa- 
ter molecules in the inner shell. Nevertheless, that g(r) 
does not rule out the possibility that the structures might 
become more flexible as the inner shell nears maximum 
capacity with lower incremental binding energies. 

To clarify these possibilities, we reduced the radius 
defining the inner sphere to R=2.68 A, for which <n> 
is close to 4 (see bottom panel of Fig. |J) and reanalysed 
the 'ab initio' molecular dynamics trajectory to extract 
the appropriate alternative moment information. Again 
using the Gibbs default model, we obtained the results 
shown in the lowest panel of Fig. |^. The inferred chemi- 
cal contribution is RTItixq « -13 kcal/mol. Using again 
the Born approximation for the outer sphere contribu- 
tion, this time with R=2.68 A, we obtain an absolute 
hydration free energy estimate of -65 kcal/mol. 

The insensitivity of these latter results to choice of in- 
ner sphere radius deserves emphasis and further discus- 
sion. From a formal point of view, the inner sphere radius 
R serves only a bookkeeping role; the left side of Eq. |l| 
should be strictly unaffected by changes in R. Neverthe- 
less, the terms on the right side of that equation are in- 
dividually affected by changes in R. Thus, we might take 
the insensitivity of the sum of those individual terms as 
an indication that the inevitable approximations are rea- 
sonably balanced. Values of R for which the sum Eq. |l| 
is insensitive are pragmatic values given the approxima- 
tions made. Fig. [7] illustrates these points and establishes 
the pragmatic value R«3.06A for the current application. 
The similarity of this value with the radius of the inner 
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shell suggested by Fig. || (3.12A) is encouraging. The 
value -68 kcal/mol is then suggested for the hydration 
free energy in the absence of any account of packing or 
van der Waals interactions. 



IV. CONCLUSIONS 

The 'ab initio' molecular dynamics simulation pre- 
dicts the most probable occupancy of the inner shell of 
Na + (aq) to be 5 and the mean occupancy to be 4.6 wa- 
ter molecules at infinite dilution, T=344 K, and a nom- 
inal water density of 1 g/cm 3 . The simulation produces 
both a satisfactory Na-0 radial distribution function and 
self-diffusion coefficient for Na + , but these satisfactory 
results required more care with thermalization and aver- 
aging time than is most common with these demanding 
calculations. Recently, this point has been separately 
emphasized in the context of 'ab initio' simulation of 
water H| . 

The complementary calculation framed in terms of 
quasi-chemical theory based on electronic structure re- 
sults for ion-water clusters, the harmonic approximation 
for cluster motion, and a dielectric continuum model 
for outer sphere contributions underestimates the inner 
shell water molecule occupancies for Na + in liquid wa- 
ter. Maximum entropy fitting of the inner shell occu- 
pancy distribution shows that the ion-water cluster re- 
sults yield a distribution significantly narrower than that 



obtained from the simulations. For this reason, naive in- 
ference of the absolute hydration free energy of Na + (aq) 
based on the cluster electronic structure results and uti- 
lizing information gleaned from the 'ab initio' molecu- 
lar dynamics was unsuccessful. The electronic structure 
calculations found minimum energy Na(fi20)7 + clusters 
only with obvious outer sphere placements of some water 
molecules though hepta-coordinate inner sphere clusters 
were observed in the 'ab initio' molecular dynamics with 
the most natural cluster definition. These results sug- 
gest that the anharmonicities and large amplitude motion 
are serious concerns, particularly for the larger clusters, 
and that the approximate theory utilized for outer sphere 
contributions must treat cluster conformations differently 
from minimum energy structures of the isolated clusters. 

Abandonment of the cluster electronic structure results 
in favor of a broader default model improved the model- 
ing of the x„ distribution on the basis of the information 
extracted from the simulation. A sequence of more ag- 
gressive fits eventually suggested the value -68 kcal/mol 
for the hydration free energy at this somewhat elevated 
temperature on the basis of the quasi-chemical perspec- 
tive of inner sphere occupancies but in the absence of any 
account of packing or van der Waals interactions. 

We acknowledge helpful discussions of many related 
issues with Gerhard Hummer and Joel Kress. This work 
was supported by the US Department of Energy under 
contract W-7405-ENG-36 and the LDRD program at Los 
Alamos. 
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FIG. 1: Structures from 'ab initio' molecular dynamics calcu- 
lations. In the top panel, the Na + ion has five (5) inner shell 
water molecule neighbors. The bottom panel shows the four- 
coordinate structure produced 70 fs later. The bonds identify 
water oxygen atoms within 3.1 A of the Na + ion. The hy- 
drogen, sodium, and oxygen atoms are shown as open, black, 
and gray circles, respectively. 
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FIG. 2: Radial distribution function gNao(r) and number n(r) 
of oxygen atoms neighboring the Na + ion. Error estimates of 
± 2a are also plotted for the radial distribution function, a 
was estimated by dividing the observed trajectory into four 
blocks of approximate duration 1.5 ps; those blocks were as- 
sumed to provide independent observations. The first mini- 
mum in the g(r) function is at r=3.12 A where g(r) falls to 
0.2. Here an average of 4.6 oxygen atoms surround the Na + 
ion. 



8 




2 4 6 8 10 12 

t(ps) 




2 4 6 8 10 12 

t(ps) 



FIG. 3: The solid line in the upper plot depicts the number 
of oxygen atoms within a radius of 3.12 A from the Na + at 
each configuration in the molecular dynamics simulation. A 
radius of 2.68 A defines the nearest oxygen neighbors in the 
lower plot. The dashed lines show the kinetic energy per 
atom during the simulation, plotted after removal of the 300 K 
thermostat at 4.3 ps. The axis on the right refers to the 
kinetic energy values. In the upper plot, an average of 4.6 
water molecules surround the Na + ion, while an average of 
4.0 water molecules surround the ion in the lower plot. 
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FIG. 4: Mean-square displacement of the Na + ion plotted 
with respect to the time interval analyzed. Analysis of the 
slope from 200-400 ps gives a diffusion constant of l.OxlO -5 
cm 2 /s. 
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FIG. 5: Free energies for Na + ion hydration in liquid water 
as a function of the number of inner shell water neighbors at 
T=344 K and ph 2 o = 1 g/cm 3 . The lowest results (open dia- 
monds) show quasi-chemical approximate values for the liq- 
uid, labelled according to the quasi-chemical interpretation. 
This graph indicates that the n=4 inner sphere structure is 
most probable under these conditions. The radius used for the 
Na + ion here is 3.1 A, though a substantial reduction of this 
value produced only a minor change in the inferred absolute 
hydration free energy; otherwise the procedure is the same 
as in previous reports [p3|, p7| . The absolute hydration free 
energy predicted here is -103 kcal/mol. The results marked 
AG (0) (filled circles) are the free energies predicted for the 
reaction Na + + n H2O in an ideal gas at p = 1 atm = p and 
T=344 K. The minimum value is at n=4. The middle graph 
(crosses) add to the ideal gas results the 'replacement' con- 
tribution reflecting the formal density of the water molecules 
-nRT ln[RTpH 2 o/p] = -n * 5.03 kcal/mol with T=344 K, 
and ph 2 o =1 g/cm 3 . 
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FIG. 6: Results for the inference of xo from 'ab initio' molec- 
ular dynamics information. The solid points represent the 
information extracted from the molecular dynamics simula- 
tion, the dotted lines are the default models, and the solid 
lines show the fit achieved by the information theory ap- 
proach. In the top panel, the hepta-occupancy was excluded, 
the probabilities were renormalized on this condition, and 
the quasi-chemical default model was used together with the 
moments (n)=4.633 and (n(n — l)/2)}=8.577. In the mid- 
dle panel, all x„ observed in the simulation were included, 
(n)=4.642 and (n(n - l)/2)}=8.624, and the ideal gas (or 
Gibbs default) model was used along with the same moments 
as above. The bottom panel shows the results using the in- 
ner sphere radius R=2.68 A and the moments (n) =4.046, and 
(n(n - l)/2})=6.393 with the Gibbs default model. 
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FIG. 7: Variation of hydration free energy contributions with 
changes in radius R defining the inner sphere. The upper 
curve is the chemical contribution obtained with an informa- 
tion theory fit using a Gibbs default. The middle curve is 
the outer sphere contribution approximated by the Born for- 
mula for a spherical ion with unit charge at its center. The 
bottom curve is the sum of the other two. The slope of the 
bottom curve is zero at R~3.06A. As discussed in the text, 
this value identifies a R-region for which the approximations 
used are pragmatically balanced. The curve takes the value 
-68 kcal/mol in that region. 



